%% Step 2

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Date: October 26, 2017
% Author: Produced for JdG,AL,IM by Christopher Evans at UPF 
%
% Program: this file reads in firm level data for values of initial import 
% shares, labor shares, and intital export shares. These values are taken
% as given (particulalry exports) in order to apply the hat matrix around a
% steady-sate.
%
% The file reads in as follows:
%   1) Import shares (gammas): for France by country (ISO coded) into an
%      array called GAMMA. Each entry into the array is a F*1 matrix. The
%      name of these files are Firm_Level_France_gamma_`ISO'.csv, and can
%      be found in the CSV/gamma_`type' folder, where `type' refers to the
%      version of imports shares imputed from the data.
%
%   2) Labor shares (alpha) and firm-level bilateral exports (pi_mnj) are
%      then imported from CSV/labor_pi.csv.
%
%   3) Domestic sales shares for French firms (pi_mmj) are then imported
%      from CSV/Firm_Level_France_pi_2.csv
%
%   4) Imports the sector list (strings), and then splits it up using 
%      unique to work out the number of firms in each sector. Then creates
%      the identifying dummy matrix used for computations later.
%
% NEW: For all code, we will change the subscripts to match those in
% master_notes3.pdf and the draft. I.e., we will change to {mn,ij} pairs,
% which represent country pair {mn} and sector pair {ij}. Global variables 
% changed and slight modification to file structure for calling up data.
% We also get rid of gamma_type variable. What previously _2.csv. files
% slightly as well as adjusting in our new directory.
%
% We also get rid of using the .(ISO{}) cell call up for the gammas and
% instead use a large matrix to save. This will help save on loops
% and speed up the procedures going forward.
%
% New: all Firm-level csv files finishing by _2 come from the "FirmModel"
% and thus treat retailing and wholesaling as NT. The corresponding files
% with no indices indeed come from the "FirmModel_WT" folder and thus have 
% firms in 51-52 into T sectors
%
% Last Updated: 5/02/2019 by IM & JDG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clearvars

global alphaT rhoT lambdaT etaT rho sigma sigmaT lambda eta varphi ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    N J FI_J france countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO ... 
    sL_n sPi_n sD_n alpha_WIOD_j alpha_WIOD_francej labour_share_PWT
 
%% Import import shares data

load MAT/WIOD.mat


% Defining France as the position in the WIOD (15th position)

france = 15;

for cty = {'AUS', 'AUT', 'BEL', 'BGR', 'BRA', 'CAN', 'CHN' ,'CYP', ...
           'CZE', 'DEU', 'DNK', 'ESP' , 'EST', 'FIN', 'FRA', 'GBR', ...
           'GRC', 'HUN', 'IDN', 'IND', 'IRL', 'ITA', 'JPN', 'KOR', ...
           'LTU' ,'LVA', 'MEX','MLT', 'NLD', 'POL', 'PRT', 'ROU','RUS', ...
           'SVK', 'SVN', 'SWE', 'TUR', 'TWN', 'USA','ROW'};
    import_file = strcat('CSV/gamma/Firm_Level_France_gamma_brn_',cty{1},'.csv');
    GAMMA.(cty{1}) = dlmread(import_file,',',1);
    clear import_file;
end

%% Import labor and export shares data

% Import the alpha, labor share scaled by WIOD

%import_file = strcat('CSV/labor_pi.csv');
import_file = strcat('CSV/labor_pi_2_brn.csv');


labor_pi = dlmread(import_file,',',1);
alpha_f = labor_pi(:,1);

% Import domestic sales shares data

import_file = strcat('CSV/Firm_level_France_pi_2_brn.csv');
pi_mnj = dlmread(import_file,',',1);

clear labor_pi

%% Import string of sectors (firms) and then create unique identifier
filename = strcat('CSV/WIOT_data/wiot12_jp_2_brn.csv');
opts = delimitedTextImportOptions("NumVariables", 2);
opts.DataLines = [2, Inf];
opts.Delimiter = ",";
opts.VariableNames = ["siren", "wiot12_jp"];
opts.VariableTypes = ["double", "string"];
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
opts = setvaropts(opts, "wiot12_jp", "WhitespaceRule", "preserve");
opts = setvaropts(opts, "wiot12_jp", "EmptyFieldRule", "auto");
wiot12jp2brn = readtable(filename, opts);
clear opts

% Work with SECTORS vector to create a dummy matrix for firms in sectors

[sector,start_f,indx]=unique(wiot12jp2brn(:,2));

%Initialises the dummy matrix that will be used for the firms in sectors

firmsecD = sparse(zeros(size(wiot12jp2brn,1),size(sector,1)));
 
for i=1:size(sector,1)
    firmsecD(indx==i,i)=1;
end

WIOD_sectors = {'AtB' , 'C' , '15t16' ,'17t18' , '19' , '20' , '21t22' ,...
    '23' , '24' , '25' ,'26', '27t28', '29' , '30t33' , '34t35' ,...
    '36t37', 'E' , 'F' , '50', '51', '52', 'H', '60', '61', '62', '63',...
    '64', '70', '71t74', 'M', 'N', 'O'}'; 

%Saving the countries to call them up easier later in loops

ISO = {'AUS', 'AUT', 'BEL', 'BGR', 'BRA', 'CAN', 'CHN' ,'CYP', ...
           'CZE', 'DEU', 'DNK', 'ESP' , 'EST', 'FIN', 'FRA', 'GBR', ...
           'GRC', 'HUN', 'IDN', 'IND', 'IRL', 'ITA', 'JPN', 'KOR', ...
           'LTU' ,'LVA', 'MEX','MLT', 'NLD', 'POL', 'PRT', 'ROU','RUS', ...
           'SVK', 'SVN', 'SWE', 'TUR', 'TWN', 'USA','ROW'};

%Adding extra row to represent end of row - makes the following loop easier

start_f(end+1)=size(firmsecD,1)+1;


% FI_J: Total number of firms (column length of firm data file)

FI_J=size(indx,1);

% firms_per_sector: Calculating the amount of firms in each sector

for j = 1:J 
   firms_per_sector(j,1) =  start_f(j+1)-start_f(j);
end

% firms_per_sector_sorted: If statement such that we have the correct order
% of sectors. Correct order is taken from WIOD. 

sector=table2cell(sector);

for i=1:J
    for j=1:J
        if isequal(WIOD_sectors{i},sector{j})
           firms_per_sector_sorted(i,1) =  firms_per_sector(j);
        else
        end
    end
end

% start_sorted: Initialise the start vector for the sorted list of sectors

start_sorted(1)=1 ; 

% start_sorted: Calculates when each sorted sector starts

for j=1:J
   start_sorted(j+1,1) =  start_sorted(j)+firms_per_sector_sorted(j);
end

% SECTORS_sorted: Ordering the sector strings in the correct order. It
%                 should follow the order of WIOD

for i=1:J
   for j= 1:J
      if isequal(WIOD_sectors{i},sector{j})
         SECTORS_sorted(start_sorted(i):start_sorted(i+1)-1,:)= ...
         wiot12jp2brn(start_f(j):start_f(j+1)-1,2);
      end
   end
end

% firmsecD_sorted: Initialising the sorted firm sector dummy

firmsecD_sorted=sparse(zeros(size(firmsecD)));

% firmsecD_sorted: Assigning value 1 for when firm is in sector j 

for j=1:J
   firmsecD_sorted(start_sorted(j):start_sorted(j+1)-1,j) = 1;
end

%Clearing old firm sector dummmy to reduce memory usage

clear firmsecD

% GAMMA_sorted. : Loops over the old GAMMAS and sorts them according to the
%                 order of the WIOD given above in WIOD_sectors 

for n = 1:N
   for i=1:J
      for j= 1:J
         if isequal(WIOD_sectors{i},sector{j})
            GAMMA_sorted.(ISO{1,n})(start_sorted(i):start_sorted(i+1)-1,:)=...
               eval(['GAMMA.',ISO{1,n},'(start_f(j):start_f(j+1)-1,:)' ]);
         else 
         end
      end
   end
end    
for n = 1:N
    GAMMA_sorted_matrix(:,(n-1)*J+1:(n-1)*J+J) = GAMMA_sorted.(ISO{1,n});
end

% pi_mnj_sorted: Initialising the pi matrix, so it can be sorted to follow
% the layout of the WIOD

pi_mnj_sorted=zeros(size(pi_mnj));

% pi_mnj_sorted: Loop that assigns the correct values to pi matrix,
% rearranging the sectors so that they are in the same order as in WIOD

for n = 1:N
   for i=1:J
      for j= 1:J
         if isequal(WIOD_sectors{i},sector{j})
            pi_mnj_sorted(start_sorted(i):start_sorted(i+1)-1,n)=...
               pi_mnj(start_f(j):start_f(j+1)-1,n) ;
         else 
         end
      end
   end
end   

% alpha_f_sorted: Making sure the firm level alphas follow the firm shares
%                 so that they are assigned to the correct firm in the 
%                 correct sector that matches up with the WIOD sector order

for i=1:J
   for j= 1:J
      if isequal(WIOD_sectors{i},sector{j})
         alpha_f_sorted(start_sorted(i):start_sorted(i+1)-1,:)=...
          alpha_f(start_f(j):start_f(j+1)-1,:);
      else 
      end
   end
end

% AlphaT = 1 : Country level labor share

if alphaT==1 
   % alpha_f_sorted: This sets the french alpha to be the same for all of
   % the firms, this is by country using the PWT
   alpha_f_sorted=ones(size(alpha_f_sorted)).*(labour_share_PWT(france));

   % Alpha_type =3 : Sector level albor share 
elseif alphaT==3
   % alpha_f_sorted: Sets the french alpha to be the same by sector for each
   %                 firm, this is from the WIOD

   alpha_f_sorted=firmsecD_sorted*alpha_WIOD_francej;
   %alpha_WIOD_j;
else
end

% If alpha_type is 2 or 4 then we use the raw Stata data for labor shares
% at the firm level. This gives us a different labor share for each firm.

%% Calculating the VA at the firm 

% Import the alpha, labor share scaled by WIOD

import_file = strcat('CSV/firm_VA_brn.csv');

% firm_VA_Turnover: importing firm level value added and turnover from a 
%                   csv file that was initially outsheeted from Stata

firm_VA_Turnover = dlmread(import_file,',',1);

% VA_f_Stata_unsorted: Saving the firm level value added 

VA_f_Stata_unsorted = firm_VA_Turnover(:,1);

% turnover_f_Stata_unsorted: Saving the turnover at the firm level 

turnover_f_Stata_unsorted = firm_VA_Turnover(:,2);


% VA_f_Stata: Sorting the VA so that firms are in the correct order

for i=1:J
   for j= 1:J
      if isequal(WIOD_sectors{i},sector{j})
         VA_f_Stata(start_sorted(i):start_sorted(i+1)-1,:)=...
          VA_f_Stata_unsorted(start_f(j):start_f(j+1)-1,:);
      else 
      end
   end
end

% turnover_f_Stata: Sorting the turnover so that firms are in the 
%                   correct order


for i=1:J
   for j= 1:J
      if isequal(WIOD_sectors{i},sector{j})
         turnover_f_Stata(start_sorted(i):start_sorted(i+1)-1,:)=...
          turnover_f_Stata_unsorted(start_f(j):start_f(j+1)-1,:);
      else 
      end
   end
end

%% Firm VA share
 
VA_f_Stata_share=VA_f_Stata./sum(VA_f_Stata);

turnover_f_Stata_share=turnover_f_Stata./sum(turnover_f_Stata);

%% Shrink large matrices with many zeros using sparse

GAMMA_sorted_matrix=sparse(GAMMA_sorted_matrix);
pi_mnj_sorted=sparse(pi_mnj_sorted);
firmsecD_sorted=sparse(firmsecD_sorted);